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1.0 INTRODUCTION AND BACKGROUND 

This report is submitted to the NASA, Marshall Space Flight Cen- 
ter's (MSFC's), Space Sciences Laboratory, in conjunction with 
the deliverable requirements of NASA Purchase Order H-08074D, 
dated August 30, 1991, entitled "Profiler Signal Detection Im- 

provements." Under the tasking provisions of this Purchase Order, 
an investigation of atmospheric signal detection algorithms was 
undertaken for the potential improvement of the predication of 
wind aloft streamline profiles. Such improved algorithm detection 
capabilities are needed in NASA's transitioning wind profiling 
radar technology for detecting atmospheric signals that are 
occasionally obscured by natural and man-made interference. 

The justification for this activity resulted from a NASA review 
of many profiler spectral samples which showed that there is a 
weak, but identifiable, atmospheric signal in the spectral power 
information of the 49.25 Mhz TYCHO Technology Inc. radar current- 
ly installed at Kennedy Space Center. This prototype system has 
an average power-aperture of 10 8 WM Z and operates with an 8 
microsecond pulse consisting of 1 microsecond code elements. 
This provides a range resolution of 150 meters over 112 range 
gates. The lowest range gate altitude of this system is located 
at 2 kilometers and the maximum range gate is located at 18.65 
kilometers. Consequently, the range gate altitudes, in meters, 
are given by the expression 2000+ ( j — 1 ) *150, for j = 1 , 2 ,..., 112 
where j equals the range gate number. 

The phased array antenna of this radar system is composed of 
coaxial -col inear elements which provide three radar beams or 
modes. The beam configuration of this system consists of two 
beams tilted 15° off the zenith on azimuths of 135° (mode 1) and 
45° (mode 2), and the third beam (mode 3) is vertical. In opera- 
tion, a pulse repetition period of 160x10 seconds is utilized 
and a cycle is completed by integrating each beam for one 
minute. For modes 1 and 2, the real time radar processor coher- 
ently integrates 320 pulses in the time domain; then, a Fast 
Fourier Transform ( FFT) is applied to this time domain data to 
obtain a 256 point (frequency domain) power spectral estimate. 
Four sets of these estimates are then averaged to provide a final 
one minute estimate. In contrast, only one power spectral esti- 
mate is formed over the entire minute for mode 3 by averaging 
1400 pulses. Consequently, modes 1 and 2 are capable of achieving 
velocity resolutions of 0.232 m/s over the ±29.6 m/s range, and 
mode 3 is capable of achieving a velocity resolution of 0.05 m/s 
over ±6.7 m/s. 

Under this current processing scheme, 256 power spectral frequen- 
cy bins are involved where the zero central frequency (bin 128) 
corresponds to the transmitted frequency. As a result, the Dop- 
pler shifted frequency bins 129-256 correspond to motions toward 
the profiler, and the frequency bins 1-127 correspond to motions 
away from the profiler. Furthermore, the TYCHO Technology, Inc. 
software supplied with the prototype system computes the radial 
range gate wind velocity component by the so-called "moment meth- 
od." In this computational method, the largest range gate power 
spectral peak is integrated over a frequency interval with inter 
val boundaries on either side of the largest peak that have power 
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amplitudes equal to the average range gate power noise level. The 
centroid frequency component of this area is taken as the Doppler 
frequency shift which corresponds to the radial velocity esti- 
mate. This computational methodology is based, on the premise that 
the power returned in this interval width is a measure of the 
degree of turbulence in the range gate volume. 

In actual practice many range gates have their largest power 
spectral peak at the central frequency, i.e., frequency bin posi- 
tion 128. For the lower range gates, this effect is often 
caused by reflections from solid objects (buildings) appearing in 
the sidelobes of the antenna pattern or local ground clutter 
effects associated with non-uniform ground heating mechanisms. 
Consequently, there are many technical issues involved. with the 
measurement of wind aloft profiles from electromagnetic energy 
scattering. The most apparent issue is that atmospheric signals 
are occasionally obscured in the presence of natural and man-made 
interference. For example, interference associated with lightning 
or aircraft detections within the scattering volume provide fre- 
quency intervals of abnormal range gate power amplitudes. In the 
case of lightning induced interference, the potential of total 
range gate power spectra contamination is theoretically predict- 
able. Aircraft detections, on the other, tend to locally contami- 
nate the power spectra altitude over time-varying frequency win- 
dows. Combining these isolated interference sources with the 
existing potential of thermal, atmospheric and cosmic noise in- 
terference provides the possibility of generating totally uncor- 
related range gate power spectra. 


2.0 PREREQUISITE ALGORITHM ASSUMPTIONS 

These implications dictate that more sophisticated data algo- 
rithms must be developed to take full advantage of this technolo- 
gy. To accommodate this situation, a preliminary investigation of 
prerequisite algorithm requirements was undertaken. During this 
investigation, the obvious question of algorithm evaluation 
criteria was raised. This prompted the following enumeration of 
prerequisite algorithm requirements and considerations: 

o Spectral data filtering should be provided by unbiased statis- 
tical estimators that account for surrounding range gate power 
fluctuations . 

o The precision of range gate spectral power measurements may 
vary because many factors have an indirect influence on these 
measurements . 

o Background noise data irregularities cannot be normalized. 

o Data filtering must be conducted to reduce or eliminate the 
magnitude of probable data error uncertainties. 


3.0 WEIGHTED MEAN OP MINIMUM VARIANCE ESTIMATOR 

By definition, a statistical estimator, u, is unbiased, if E(u) = 
0. This property merely states that the random variable u pos- 


2 



sesses a distribution whose mean is the parameter 8 being esti- 
mated. To place this concept in the context of the current prob- 
lem, suppose that ^ denotes the power amplitude of the j 
range gate and i^*"* frequency, for a selected beam. Here, i — 
1 , 2 ,..., 256 and j = 1,2,..., 112. Now restrict j to a set of 2n+l 
range gate amplitudes, by taking n range gates on either side of 
the k*'® range gate of interest. By letting j = k-n,k-n-l, . . . , 
k, . . . ,k+n-l,k+n one obtains 2n+l power amplitude range gate ob- 
servations for the i th range gate frequency. If each of these 
observations are assumed to have the same expectation, one can 
assume that subsequent range gate velocities are linearly corre- 
lated. This condition is imposed by considering the following 
linear combination of 2n+l random variables: 


k+n 


Since E(u.s) must equal 0 for u A to be unbiased, it is customary 
to impose a measure of the variability of about 9. This can be 
accommodated by forcing the condition of minimum variance. To 
Introduce this property of unbiasedness and minimum variance, let 

k+n k+n k+n 

E(Ui) = E ( Sa-jX.-j) = Ea j E(X i j)= 9 2 a 

j=k-n J j=k-n 3 =k-n 

This formulation introduces the constraint that 

k+n 

Ea.: = 1 
j=k-n 


since E (U-: ) must equal 0 for to be unbiased. Furthermore, the 
variance of u^ is given by the following expression: 


V(Ui) 


k+n 
V( Ea 


l x i i ) 
j=k-n X '^ 


k+n 


. sa j V(X i,j } 

j=k-n 


k+n „ 

= Ea^ 2 a 
j=k-n 



To determine the values of the a.: coefficient weights which make 
the variance of u^ a minimum, one may write the sum of these 
weights in the following form: 


a k+n + 


k+n-1 

Ea-: 

j=k-n 


k+n 

Ea-s = 1; 
j=k-n 
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so that 


a k+n 


k+n-1 

j=k-n 


Thus, the variance of can be written in terms of the first 
2n+l weights by the expression. 


V(Ui) 


k+n_ J 2 

= 2a + 

j=k-n J 


a k+n a k+n 


k+n "i 2 

= Za-^a-s'* + 

j=k-n J 


k+n-1 
[i- z ai ] 
j=k-n 


2 


a k+n 


2 


k+n-1 
= Za.t 2 ff 
j=k-n 


j 


2 


+ 


k+n-1 k+n-1 

1 - 2 Za-: + ( Za-j ) 2 ]ff k+n z . 
j=k-n j=k-n 


To find the values of a.: which make this variance a minimum, one 
must differentiate with respect to aj and equate to zero: 


dV(fL i ) 


da. 


2 [ aj<Jj‘ 


k+n-1 

- ( 1- z aj )a k+n 2 ] 
j=k-n 


- 2ajUj 2 - 2a k+n a k+n 


= 0. 


Upon solving this last equation for a.: and letting j take on the 
values k-n,k-n+l, . . . ,k, . . ., k+n-1 one obtains the expression 


0 _ 2 
a k+n a k+n 



However, this relationship is obviously true for j = k+n also. 
Consequently, one can sum over j = k— n,k— n— 1, . . . ,k+n to obtain 


k+n 


k+n 


2 a j - a k + n a k+n 2 

j=k-n D=k~n 
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Since the sum of these aj coefficient weights equals 1, 


k+n 

<Jj 2 S Wj 

J j=k-n 


for j = k-n,k-n+l, . . . , k, . . . ,k+n-l,k+n, where Wj = l/Uj 2 . Upon 
substituting these coefficient weights into the above variance 
equation, 


V(tii) = 


k+n 
2 W-j 
j=k-n 


equals the minimum variance. Here one must realize that o j 
equals the j® 1 range gate power amplitude sample variance. 


m 


. s < x i,j “ Xm > 
i=l ' J 


m-1 


where m = 256 and equals the arithmetic mean of the Xj^ j power 
amplitude observations, i.e., 



m 


In su mm ary, this unbiased weighted mean of minimum variance for- 
mulation provides the capability to statistically construct a 
mean range gate power spectrum which accounts for adjacent range 
gate data fluctuations. The compromise associated with this ap- 
proach resides in the assumption that the wind velocity is as- 
sumed to be linearly correlated over subsequent range gate vol- 
umes. However, the main disadvantage associated with this formu- 
lation is that the magnitude of persistent non-canceling range 
gate power amplitude fluctuations and irregularities are not 
systematically removed or filtered out. That is, this formulation 
does not reduce or eliminate non-canceling range gate probable 
data error uncertainties. Consequently, additional data filtering 
or smoothing is required. 
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4.0 DATA FILTER SELECTION 


To help understand the fundamental data filtering requirement, a 
literature search was conducted of approximately 43 different 
sources, which revealed two general points of interest: 


o Moving average data filters contain empirical elements that do 
not appear to be logically necessary. 

o The problem of data filtering belongs to the mathematical theo- 
ry of probability. 


4.1 MOVING AVERAGE FILTERS 

The utilization of moving average data filters, such as simple 
moving average or weighted exponential moving average, contain 
empirical elements that do not appear to be logically necessary. 
Specifically, each of these filters are based on the assumption 
that the true signal can be expressed over a finite time window 
by a preselected functional curve type. For instance, the simple 
n point moving average filter gives the average of the n most 
recent observations and is equivalent to fitting a first degree 
polynomial (line) to the observations evaluated at the observa- 
tion window midpoint. If the process is constant and n is large, 
the simple moving average estimate will be stable. However, if 
the process is changing, a small value of n is needed for rapid 
response. Furthermore, each of these moving average filters have 
a graduated weighting scheme over their finite data smoothing 
window. Apart from mere convenience, these empirical elements are 
theoretically questionable. 


4.2 PROBABILISTIC FILTER IMPLICATIONS 

It is interesting to note that the arbitrary selection or need 
for these empirical elements can be removed if one views the 
filtering process in terms of mathematical probability theory. 
One such probability approach is to assume that the data observa- 
tion irregularities can be viewed in terms of the hypothesis that 
data observation uncertainties are distributed in accordance with 
the normal probability law of errors. In theory, this statistical 
graduation methodology has the capability to utilize all of the 
available observations to obtain the "most probable" set of ob- 
servations, in comparison to the selection of an arbitrary data 
window smoothing domain. In situations where only erroneous iso- 
lated peaks are present in the data, this approach should provide 
extremely reliable signal reconstruction estimates. This is a 
direct result of the fact that these isolated peaks are a low 
probability event. By comparison, the situation of having random 
errors in the data observations (i.e., errors of zero mean with 
serially uncorrelated root-mean-square values) should also pro- 
vide for extremely reliable signal reconstruction. This perform- 
ance prediction is a direct result of the inductive antecedent 
probability that if the observations had been more accurate, the 
unfiltered data observations would have been smooth. 
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5.0 NORMAL PROBABILITY LAW FILTER DEVELOPMENT 

Under this proposed filtering approach, the basic provisions of 
the Fundamental Theorem of Inductive Probability apply . Unfortu- 
nately, this theorem imposes the embedded assumption that the 
measurement precision is the same for all data observations. To 
understand the physical implications of this assumption, suppose 
one is concerned with n equidistant data observations, u n , that 
are affected with uncertainties or irregularities due to acciden- 
tal errors of observation. In other words, when the u n data is 
plotted as a function of n, the points do not lie on a smooth 
curve, although there is a strong antecedent probability that if 
the observations had been more accurate the curve would be 
smooth. Here, the term "smooth" may be given a more precise defi- 
nition by interpreting it to mean that the third order forward 
data differences, S Uj = Uj— 3Uj_ 3+3Uj_2~ , must be small. 

To mathematically express this concept of smoothness, the nota- 
tion of forward differences is used. This notation is based on 
the operator 5, which expresses the differences of equally spaced 
data observations. In this notation, the coordinates of a point 
are given by (t k ,u k ) where the t k abscissa are equally spaced 
values with h delta spacing (i.e.,t k+ n = t^+(k-l)h) and the ordi- 
nates, u k , are expressed in the following functional notation, 
u(t k ) = u k . In this notation, the first forward difference of a 
tabular point u(t k ) is defined as 

6u(t k ) = u(t k +h) - u(t k ) or Syl = u k+1 - u k . 

As a result, forward differences of various orders are given by 
the following general expressions: 

^k 1 = u k+l “ u k' 

5 k 2 = u k+2 “ 2u k+l + u k' 

5 k 3 = u k+3 “ 3u k+2 + 3u k+l ” u k' 

<S k 4 = u k+4 - 4u k+3 + 6u k+2 “ 4u k+l + u k' 


n n! 

6 k n = ZC-lp [ 

j=0 ji(n-j)! 


u k+n- j * 


Now consider the hypothesis: 

The true value, which should have been obtained by the 
observation u^, lies between u. and u^ +*, where a 
is a small constant (e.g.,one unit in the last decimal 
place used in the data observation measurements) , so 
that the true value which should have been obtained by 
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the measurement for u 2 , lies between u 2 ^ and u 2 ?"+a, 
and so on, up to u n ; which lies between u n x and u n x +A. 

This shall be called "Hypothesis H." However, before the observa- 
tions have been made, one has no probability knowledge of the 
situation, other than the above hypothesis and the definition of 
smoothness. Nevertheless it is possible to formulate the degree 
of smoothness of the sequence of Uj measurements, by taking the 
sum of squares of the third differences 

S = (u 4 1 -3u 3 1 +3u 2 1 -u 1 1 ) 2 + (u 5 1 -3u 4 1 +3u 3 1 -u 2 1 ) 2 + ... + 
( u n 1 - 3u n-l 1+3u n-2 1 “ u n-3 1 ) 2 

which will be large for sequences of rough data. Therefore, S may 
be called the measure of roughness of the sequence. In theory, 
this measure can be extended to the case where the observations 
are not taken at equidistant values of time by taking, instead of 
S, the sum of the squares of the third divided differences of the 
graduated values. 

By analogy to the normal probability law of errors, one may sup- 
pose that the a priori probability of Hypothesis H is 

Cexp(-B 2S )A n (5.0-1) 

where C and B are constants. Next consider the a priori probabil- 
ity that the observations will be u 1 ,u 2 ,u 3 , . . . ,u n , under the as- 
sumption that Hypothesis H is true. 

Since the true value of the first observed quantity is u-,^ 1 (under 
this hypothesis) , the probability that a value between u^ and 
u-^+a will actually be observed is 

hj/. J ir expf-h-j 2 ^ a, 

under the postulated normal probability law of errors, where h^ 
is a constant which measures the precision with which this obser- 
vation can be made and = (u^-u^ ) . Similarly, the probability 
that a value between u 2 and u 2 +* will actually be obtained for 
the second observed value is 


h 2 /7ir exp(-h 2 2Qt )A 

where a = (u.,-^ 1 ) 2 and h 2 is the measure of precision of this 
second observation. Thus, under the assumption that Hypothesis H 
is true, the a priori probability that the first observed quanti- 
ty will lie between u-, and u^*, the second observed quantity be- 
tween u 2 and u 2 +a, and so on, is 

[h 1 h 2 h 3 . . .h n /(7(r) n ] exp(-n) * n (5.0-2) 

where Cl = h-j^ (u 3 — u^ 3 ) 2 + h 2 2 (u 2 — u 2 ^) 2 + ... + h n (u n — u n ) 

The sums S and n enable one to numerically express the smoothness 
of the graduated values, as well as the fidelity of the graduated 
to ungraduated values . 
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To accommodate this formulation, the Fundamental Theorem of In- 
ductive Probability must be utilized. That is, suppose that a 
certain observed phenomenon may be accounted for by any one of a 
certain number of hypotheses. Then suppose that the probability 
of the s th hypothesis (based on a priori information before the 
phenomenon is observed) is p s . Under this situation, the proba- 
bility of the observed phenomenon, based on the assumed truth of 
the s th hypothesis, is P g . Consequently, the probability of the 
s'"* 1 hypothesis is 

Ps^s/^ (Ps j ^*s j ) 

where S denotes the summation over all the hypotheses j = 
1,2,3,... . It follows from this expression that before the phe- 
nomenon was observed the most probable hypothesis was that for 
which p g was the greatest. However, the most probable hypothesis 
after tne phenomenon has been observed is that for which the 
product P S P S is the greatest. Applying this condition to the case 
under consideration, equations (5.0-1) and (5.0-2) may be com- 
bined to obtain 

[Ch 1 h 2 h 3 . . .h n /(y r ) n l exp(-B 2S_n ) A 2n . 

o ^ — n . . 

Since this expression is a maximum when is a minimum, the 
most probable set of values, u-, ,u 7 ,u-» ,...,u_ , are obtained 
when B 2S-n is a minimum. 


5.1 ANALYTICAL MINIMUM FORMULATION 

Writing down the ordinary difference equation conditions for a 
minimum, one obtains the equations: 


h 

h 

h 

h 


1 

2 

3 

4 



^l 1 

- 

BVu^ 

h 2 2u 2 1 

+ 

3B 2 <S 3 u 2 1 - B 2 5 3 u 2 3 “ 

h 3 2 u 3 1 

- 

3B 2 «S 3 u 1 1 + 3B 2 i 3 u 2 3 ’ - B 2 6 3 U 3 1 

h 4 2 u 4 1 

+ 

B 2 <S 3 u 1 1 - 3B 2 6 3 u 2 1 + 3B 2 <5 3 u 3 1 - B 2 <S 3 u 4 1 


h n 2u n - h nV 


B 2 «S 3 u 


n-3 


By making the simplifying assumption that the measure of preci- 
sion is the same for all data observations, i.e. ( h^ = h 2 = h 3 
= ... =h n ) , one imposes a constant measure of precision across 
all observations. If this is not the case, one can graduate some 
function of u, such as ln(u) , instead of u, choosing this func- 
tion so that its measure of precision has nearly the same value 
for all values of the argument. However, under the equal preci- 
sion assumption, one may substitute hj = cB 2 for j = l,2,3,...,n 
into each of the above equations. 
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Here, e denotes the index of precision generally associated with 
the normal probability density function. 


6 o 2 

f (t) « exp(-e^) . 

J IT 

(See Section 5.5, NORMAL PROBABILITY INDEX OF PRECISION FORMULA- 
TION, for a detailed discussion of this function.) 

Under the above noted substitution the equations have the form. 


eu 3 = 

eUf 1 

- 

S 3 ^ 1 

eu 2 = 

6U2 1 

+ 

36 3 u 1 1 - 5 3 U 2 1 

eu 3 = 

eu 3 1 

- 

3d 3 u 1 1 + 35 3 U 2 1 “ ^u-j 1 

eu 4 = 

6U4 1 

+ 

fi 3 U 1 1 - 36 3 u 2 1 + S^u-j 1 - 5 3 u 4 1 


eu n " 6U n X + 63u n-3 1 * 


Now all of these equations, except the first three and last 
three, are of the form 


€Uj = euj + ^Uj.j 1 . 

This can be directly verified analytically from the forward dif- 
ference equations provided earlier. If one introduces a quantity 
u 0 \ such that S^Uq = 0, the third equation, 

eu 3 = eu 3 ^ - 5 6 Uq 1 , 


is of this same form. Similarly, the first two and last three 
equations can be brought to this form, by introducing the new 
quantities / u n+l ' u n+2 • an< * u n+3 ' suc ^ that <S u_ 3 = 

0, 5 3 u_ 2 1 = 0, «S 3 u n _ 2 1 = 0, «S 3 u n _ 1 1 = 0, and S 3 ^ 1 = 0. Thus the 
graduated values Uj 1 satisfy the linear difference equation 

euj 1 - 5 6 Uj_ 3 1 = euj , (5.1-2) 

being in fact the particular solution of this equation which 
satisfies the six terminal conditions 



6 3 u_ 2 1 = 0, S 3 u _ 1 1 = 0, S 3 m_ 2 1 = 0, $ 3 n _ 2 - 0, 5 3 u n _ 1 1 = 0, and 

6 3 u n X = 0, 

which provides the following relationships 

6 4 u _ 2 1 = 0 , 5 4 u _ 1 1 = 0 , <S 5 u _ 2 1 = 0 , 5 4 u n _ 2 1 = °» ^ 4 u n -l 1 = °» and 
S 5 u n. 2 1 - 0 (5.1-3) 


5.2 CONSERVATION THEOREM OF MOMENTS 
By the summation of equation (5.1-2) 

e (u 1 1 +u 2 1 +. • .+UH 1 ) -€ (u 1 +u 2 +. . * +u n ) = 6 6 u_ 2 1 +'5 6 u_ 1 1 +. . ,+ 6 6 u n _ 3 1 

- x5„ 1_*5„ 1 

~ 5 u n -2 5 u -2 

= 0 , 

as a result of the 3 rd and 6 th particular solution values of 
equation (5.1-3). 

Therefore, 

u i 1 +u 2 1+ . * * +u n 1 = u i +u 2 + *** +u n * 

Moreover, by equation (5.1-2) 


6 SjU-j 1 - 

j=l 


or u^ 1 + 2 u 2 


e Eju.: = ^(k+S )^^ 1 = 0 

j=l J k =-2 

1 +...+nu n 1 = u 1 +2u 2 +. . .+nu n . 


At this point it should be noted that equation (5.1-2) can be 
written in the form 

€j 2 (Uj 1 -Uj) = j 2 5 6 Uj_ 3 1 

for j = 1 , 2 , 3 , ...,n. Upon 'expanding over j and adding, these 
difference equations, 


n _ . n , n_J r 

e Zj 2 u-: 1 - e Zj 2 Ui = S(k+3) 2 «S 6 u k i 

j = l J j=l J k =-2 

= <S 6 u_ 2 1 +2 2 «S 6 u_ 1 1 + +n 2 5 6 u n _ 3 1 



= n 2 <S 5 u n _ 2 1 -(2n-l)6 4 u n _2 1 +25 3 u n _ 2 1 -6 3 u_ 2 1 - 


= 0 . 


Therefore, Ui ^+ 2 2 u 2 ^+ . . . +n 2 u n ^ — 11^+2^115+. . ,+n^u n , which shows 
that the moments of order 0,1, and 2 are the same for the filtered 
and unfiltered data. In other words, the filtered and unfiltered 
graphs have equal area, center of gravity, and moments of inertia 
about any line parallel to the observation ordinates. 


5.3 GRADUATED DIFFERENCE EQUATION SOLUTION 
In this section, the solution of the difference equation 

euj 1 - S 6 u j.g 1 = euj, 

subject to the six terminal conditions outlined in Section 5.1 
will be investigated. To start this discussion, assume that the 
normal probability index of precision, e, is given, so that a 
general solution can be obtained in terms of symbolic operator 
notation. Here, it is assumed that e is a positive non-zero con- 
stant. Since 6 = E-l in operator notation, this difference equa- 
tion can be written as: 

[ { (E-l) 6 -eE 3 }/E 3 ]Uj 1 = -euj 

or Uj 1 = -e[E 3 /{ (E-l) 6 -eE 3 }]Uj ( 5 . 3 - 1 ) 

In this form, considerations of symmetry lead one to expect that 
each u 1 will be given as a linear function of the u's, with coef- 
ficients symmetrically placed about the central term. This sug- 
gests expanding th© right-hand side in powers of both E and E . 
Such an expansion is readily seen to be possible, because the 
complex equation (z— 1) — ez — 0 has its six roots in reciprocal 
pairs, so that three are less in modulus than unity. Hence, the 
equation [ ( z- 1 ) °-ez 3 ] _1 may be expanded as a Laurent series, 
convergent for z = 1, and the coefficients may be found by the 
usual theory. The fact that these coefficients are symmetrical 
follows from the fact that the operator is left unaltered by the 
substitution of E -3- for E. 

Utilization of the theory of functions of a complex variable may 
be avoided by having recourse to the theory of partial fractions. 
That is, if the operator in question is resolved into six partial 
operators by this theory, three of these operators may be expand- 
ed in powers of E and the other three in powers of E . As a 
result, the addition of these operators provide the complete 
expansion. 

Either of these methods leads to the same graduation formula, 
namely : 

Uj 1 = k 0 Uj + ki( u j+i +u j-i) + k 2 (Uj + 2 +Uj_ 2 ) + ..., ( 5 . 3 - 2 ) 


1 O 


where 


< Sl ) n+2 

n (Si-S 2 > (Si-Sj) (s^CSjT 1 ) (s 1 -(s 2 )" 1 ) (Sj-fSj)' 1 ) 

Here Sn ,s 2 ,S 3 are the three roots less in modulus than unity of 
the equation (z-lj^-ez"* = 0, and E denotes interchange of 
s l' s 2' s 3' followed by summation. One of the roots is evidently 
real; hence, these roots may be written as r^, r 2 exp(i0) , and 
r 2 exp(-i0) , so that the above equation becomes; 


k j = " 


r 1 2 -2r 1 r 2 Cos[9]+r2 2 


. 5-s 3 . {r 1 ^~ 1 H 


Q j 

Sin (0) 


} 


where 


Qj = r 2 ^” 2 . [r 2 Sin( (j-2)0)-r 1 Sin((j-l)0) ] 
for ] = 1,2,3, 


This expression provides the necessary provisions needed to com- 
pute the k-s coefficients needed in the graduation formula of 
equation (5\ 3-2) for representative values of the normal proba- 
bility index of precision. 


5.4 EXTENDED AUXILIARY DATA POINT PROVISIONS 


An inherent defect is present in this formulation unless some 
means can be developed to graduate the complete set of data ob- 
servations. To investigate this possibility, consider attaching 
auxiliary data points to each end of the observation data set 
u. , u 2 , u 3 , . . . , u R . The existence of this possibility is implied by 
the graduation formula of equation (5.3-2). To accommodate this 
addition of auxiliary data points beyond the domain of the origi- 
nal data observations, the requirement can be imposed that the 
third differences involving these extended points must be zero. 
Insight into this formulation can be obtained by first introduc- 
ing the auxiliary points, u n+; , ' u n+2 ,u_ +3 , . . . and u 0 ,u_n ,u_ 3 , . . . . 
By differencing the added tnird differences three rimes, one 
obtains 


<S^u n _ 2 ^ = 0, ^ u n-l^ = 0 , ... , and 6^u_ 3 ^ — 0, <S^u_ 4 ^ — 0 ... . 

Referring to the difference equation (5.1—2), and imposing the 
condition that = u^ and u^ = u^ for i = 0,-1, -2,... and k = 
n+l,n+2,n+3, . . . forces the condition of having extended^raduated 
and ungraduated values coincide. It also follows that 6 u n+ ^ = 0, 
5 3 u n+2 =0, ..., and 5 3 u_ 3 = 0, 6 3 u_ 4 =0, ... , as a result of 
this condition. Consequently, equation (5.1—2) along with its six 
terminal conditions combined with the above auxiliary conditions 
of having the indicated third differences equal to zero, suffices 
to determine the extended auxiliary data requirements. 
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If the original number of data observations is not too small 
(i.e., if it exceeds twenty- f ive) , the extended data may be found 
as accurately as required by the following method. From the con- 
dition that the graduated and ungraduated auxiliary data coincide 
together with equation (5.3-2) provides 


u 


j 


eE 


3 


(E-1) 6 -€E 3 


u 


j 


(E-l) 6 

or u-s = 0, for j = n+1, n+2, ... . 

(E-l) 6 -eE 3 J 

This last equation may be written as, 

(E-l) 6 (E-l) 3 

. — I.Uj = 0. (5.4-1) 

(E-s 1 ) (E-s 2 ) (E-s 3 ) (E-ts-L)" 1 ) (E-(s 2 ) 1 ) (E-(s 3 ) - 1 -) 

To show that u_ +1 , u n+2 , ... are given by the equation formed by 

retaining only the first operating factor on the left side of the 
above equation, namely; 


(E-l)-* 

(E-s-l) (E-s 2 ) (E-s 3 ) 


Uj = 0 


(5.4-2) 


one should note that this equation gives Uj 

u-j_ 1 , u^ — 2 / / since the operator can be 

ing powers of E, in the form 


uniquely in terms of 
expanded in descend- 


1 — q^E 3 — q 2 E 3 — g 3 E — ... , 


so that 


Uj = qiUj_ 3 + q 2 u j-2 + ^3 u j-3 + ••• ' 

for j = n+1, n+2, n+3, ... . From this equation u n+1 , u n+2 , ... 

are determined in succession, provided that the terms containing 
the undetermined extended data points Uq, u_ 3 , u - 2 ' ••• are ne 9“ 

ligible. These values of Uj will satisfy equation (5.4-1) for the 
operator _ 

(E-l) 3 

(E-( Sl ) _1 ) (E-fSj)*" 1 ) (E-(s 3 )‘ 1 ) 

which is expandable in ascending powers of E. Now equation (5.4- 
1) can be deduced from equation (5.4-2), which holds for j = n+1, 
n+ 2 / 3 ^ • • • • 

Lastly, the conditions <S 3 u n+1 = 0, <S 3 u n+2 = 0, ... , and <5 3 u_ 3 = 

0, 6 3 u_ 4 =0, ... will be satisfied since equation (5.4-2) can be 



written as 


(E-l) 3 


( E-s ± ) (E-s 2 ) (E-s 3 ) 


(E-s^ (E-s 2 ) (E-s 3 )Uj = 0. 


Hence (E-1) 3 u^ = 5 3 u-s =0, for j = n+1, n+2, n+3, ... . In addi- 
tion, the values of gj for j = 1,2,3,... are found to be 


q j = ’ 


6 • • 
6-j 


1 

— t v j+3 - 3v j+2 + 3v j+l - V jJ 


where D = r^ 2 - 2r^r 2 Cos[9] + r 2 and 


Vj = rj 


. r i^” 2 {r 2 Sin( ( j -2 ) 6) + r 1 Sin( (j-l)9) } 

Sin(9) 


where r^, r^exp(i9) and r 2 exp(-i9) are the real and complex roots 
of the equation 


( Z— 1) 6 - €Z 3 = 0. 

a f Vi 

To compute these roots, expand (z— 1) using the r term expan- 
sion formula of the Binomial Theorem, 

— a n_r+1 b r ” 1 of (a+b) n 

(r-1) ! (n-r+1) ! 

where a=z, b — — 1 an n = 6 for r = 1,2, 3, 4, 5, 6, 7. That is, 

r = 1 — z 6 (-1)° = z 6 

016 ! 


r = 2 


6! 5 1 5 

Z 5 (-1) 1 = -6z b 

1151 


z 4 (-1) 2 = 15z 4 

214 ! 


15 


u . 

r = 4 z 3 (-1) 3 = -2 Oz 3 

3 ! 3 ! 


6 ! 

r = 5 z 2 (-1) 4 = 25z 2 

4 ! 2 ! 


y • 

r = 6 z 1 (-1) 5 = -6z 

5111 


6 ! 

r = 7 z° (-1) 6 = l. 

6 ! 0 ! 


As a result, 

(z-1) 6 - ez 3 = z 6 - 6z 5 + 15z 4 - ( e+20) z 3 + 15z 2 - 6z + 1 = 0 

where r-, equals the smallest real root of this equation and r 2 
equals the modulus (i.e., r 2 = ./(c^ 2 +c 2 2 ) ) of the complex conju- 
gate root pair, z = c-^ ± ic 2 , less than unity and 

9 = Tan -1 ( 02 / 0 -^) . 

Since the successive values of q_j decrease quite rapidly the 
number of auxiliary data points needed in practice will depend 
upon the degree of accuracy required and the rapidity with which 
the successive graduating coefficients diminish. In theory, only 
three auxiliary data points at each end need be calculated. The 
rest can be found by virtue of the constraints that (S 3 u^ = 0 and 
S 3 u k = 0 for j = n+1, n+2, n+3, ... and k = -3, -4, -5, ... , by 
forming a difference table with zero third differences. Notice, 
exactly the same process may be applied t-o both ends. However, 
the far end may be handled by reversing the order of the original 
data observations. 


5.5 NORMAL PROBABILITY INDEX OF PRECISION FORMULATION 

As noted earlier, representative values of the normal probability 
index of precision, e, may be used to compute the graduated sig- 
nal reconstruction coefficients, k^, as well as, the auxiliary 
data point coefficients, q k . The underlying consideration is the 
determination of which values of e are likely to be of practical 
application. To consider this dilemma, it is necessary to find a 
means of assigning, in advance, an appropriate value of e. Sta- 
tistically, this value should be based on the probable error of 
the sampled data observations, because many unknown factors have 
an indirect influence on these values. In addition, such a sta- 
tistical approach has the advantage of preserving the probabilism 
of this formulation. 
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Under this approach, one may utilize the analytical form of the 
normal probability *index of precision integral, discussed by 
Worthing and Geffner : 


f (t) dt = 


f (t)dt 


-Z 


1) The total number of deviations, t = u - n u , over the 
entire range of possible values of t must be equal 
to the total number of observations, n, i.e.. 


00 


n = n 


f (t) dt or 


f (t) dt = 1. 


—00 


where n u equals the average value of the original u 
observations 


2) 

3) 


1 n 


n u = 


Z u- 

n j=l ‘ 


The most probable value of t = 0, i.e., 


u = n u . 


When u = n u or t = 0, f(t) is a maximum. 


the value 


To determine the analytical form of f(t), assume that the values 
of t can be divided into At intervals so that only one value of t 
occurs in each interval: 


f(t x ) At 
f (t 2 ) At 


probability of the deviation t-^; 
probability of the deviation t 2 ; 


f (t n ) At 


probability of the deviation t n . 


* Worthing, A.G. and Geffner, J.; "Treatment of Experimental 
Data," John Wiley and Sons, Inc., New York (1943), p.155. This 
integral is also listed as No. 492 in Pierce, B.O.; "A Short 
Table of Integrals," Ginn and Co. (1929). 
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Therefore, the probability of this combination of deviations is 
given by the expression, 


n 

W = (*t) n E f (t-j) . 

j-1 


Upon taking the logarithm of both sides of this equation: 

ln(W) = Zln(f (tj ) ) + nln(At) (5.5-1) 

From restriction 2 and the definition of n u it follows that 

n n n 

St., = Z(Uj - n u ) = ZUj - nn u = 0. 

j-1 j=l j=l 

Given this result together with restriction 3), it follows that 

din (W) 

= o 

dh u 

so that the probability W depends indirectly on n u through the 
t-'s. By utilizing indirect differentiation (i.e., the Chain 
Rille) on equation (5.5-1): 

d(lnW) d(lnW) dt-L d(lnW) dt 2 d(lnW) dt n 

= . + . + . . . + “ - 0 • 

dn u dt-L dn u dt 2 dn u dt n dn u 

Since the last term of equation (5.5-1) is a constant, its deriv- 
ative is zero, so that each term of this equation involves only 
one t. This implies that, 


d(lnW) n 

= Eg (t-s ) = 0 

dn u j-1 


(5.5-2) 


where a typical term is defined by the expression 

df(tj) 


dtj 


g(tj) = 


f(tj) 


(5.5-3) 


In arriving at this last result, the fact that 


dtj d (Uj -n u ) 
dn u dn u 

for j = 1,2,3, ...,n was utilized. In line with the standard tech- 
niques of solving differential equations, each of these g(tj) 
terms can be expressed by an infinite series. 
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g ( t j ) = bp + b^t j + b 2 tj ^ ^ ••• ••••(5.5 4) 

These constant coefficients may be determined from equations 
(5.5-2) and (5.5-4) by summing the respective g(tj) series for j 
= 1,2,3, ...,n and setting the resulting sum equal J to zero in ac- 
cordance with the constraint of equation (5.5-2). These opera- 
tions provide, 

n n n n 

Zg(ti) = nb n + b St^ + b 2 Zt^ + b 3 Zt-: J + ... =0. 
j-1 3 ° 3 = 1 3 j=l D j-1 3 

In addition, the above equation can be satisfied if each of the 
individual terms are equal to zero. This may be accomplished by 
having all of the b's except b^ equal to zero. This condition is 
imposed because 


n 

St-: = 0. 

j=l 

By utilizing equations (5.5-3) and (5.5-4), in conjunction with 
the above bj coefficient constraint 

df(t i ) df(t-s) 

— = bit.: f (t^ ) or = b x tjdtj . 

dtj 3 J f(tj) 

Upon integration, 


df (t) 
f(t) 


= In ( f (t) ) + C, = b. 


b, t‘ 


tdt = 


+ c 2 . 


b x t 


and In ( f ( t ) ) = + (C 2 - C x ) or f(t) = C exp(b x t /2) 


where C is the constant of integration. Since f(t) is known to 
have a maximum at t = 0 , and is expected to approach zero for 
very large positive or negative values of t, b 3 must be negative. 
Thus, it is appropriate to replace b-^/2 by 


€ 


2 


b 


1 


2 


which yields f(t) = C exp(-e 2 t 2 ). The constant of integration, C, 
may be evaluated by substituting f(t) into restriction 1) and 
performing the indicated integration, i.e., 


00 


oo 


f (t) dt = 2 


f (t) dt = C 


exp(-e 2 t 2 )dt = = 1 


c7(r 


— 00 


-oo 
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Therefore, the constant of integration is given by 


C = and f(t) = exp(-e 2 t 2 ) 

y it j ir 

As a result, a large value of e provides a high, narrow curve for 
f(t), which corresponds to observations with a high degree of 
precision. Conversely, a small value of e yields a low, flat 
curve which corresponds to observations of low precision. 


5.6 PROBABLE OBSERVATION ERROR FORMULATION 


To formulate the probable observation error, one must establish 
the probability of obtaining an observation with deviation, t, in 
the range -P e to P e . This probability is found by integrating 
f(t) from -P to P e and setting this integral equal to ^ . Howev- 
er f(t) is a symmetric function about t = 0 so that this integral 
is twice the integral from 0 to P e . Hence, the probable error, 
P e , of an observation is defined as, 


r 

r e 26 r 

f ( t) dt = 2 

f ( t ) dt = — 

J 

J (T . 

(D 

O 

O 


exp(-e 2 t 2 )dt = 


In other words, P e is that magnitude of deviation whose probabil- 
ity of being exceeded is one-half. This consideration gives: 


P 


e 


0.476900036 


e 


The geometrical interpretation of P e is that the area under f(t) 
between the limits -P e and P e is equal to h the total area under 
f(t). 

Another interesting relationship is that a = \/ej2, where a is 
the standard deviation of the original observations. To derive 
this relationship, rewrite the standard equation for variance in 
the form 


a 


2 


00 

e f 


J r J 

—00 


t 2 exp (-e 2 t 2 ) dt 


00 

2 e f 


J IT J 
0 


t 2 exp (-e 2 t 2 ) dt 
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00 


ej ir J 


texp(-e 2 t 2 ) (-2e 2 t) dt 


and integrate by parts using 
exp(-e 2 t 2 ) to obtain 


U = t, dV = - 2 te^exp(-e 2 t 2 ) and V = 


a 


2 


1 


ey IT 


00 

* 

UdV 

0 


| exp(-e 2 t 2 )dt) . 

J 

0 

Clearly, the first term on the right vanishes at both limits and 
the second term is equal to J (r/2e, so that 

1 1 0.707106781 

a 2 = or a = ~ • 

2e 2 ej2 £ 


oo 


e7r 


{ [texp(-e 2 t 2 ) ] 


5.7 REPRESENTATIVE GRADUATING COEFFICIENTS 

By means of the above methods, the graduating coefficients, k n , 
and auxiliary data coefficients, q p have been evaluated for a 
number of representative values of the normal probability index 
of precision, e. Tables of these values to four decimal places 

are given below: 


k n Coefficients 

Graduation Formula: 

Uj 1 = k 0 Uj + k i( u j+l +u j-l) + k 2 (Uj + 2 +Uj_ 2 ) + ••• • 

n e=0.01 €=0.02 e=0.05 e=0 . 10 €=0.25 €=1.00 


0 

1 

2 

3 

4 

5 

6 

7 

8 
9 


0.1570 
0.1482 
0 . 1254 
0.0948 
0 . 0628 
0.0341 
0.0117 
-0.0035 
-0 . 0119 
-0.0148 


0.1769 
0.1644 
0.1329 
0.0928 
0.0536 
0.0216 
-0 . 0004 
-0.0125 
-0.0165 
-0.0151 


0.2076 

0.1873 

0.1397 

0.0842 

0.0359 

0.0027 

-0.0145 

-0.0191 

-0.0161 

-0.0099 


0.2347 

0.2056 

0.1412 

0.0721 

0.0191 

- 0.0110 

- 0.0211 

-0.0186 

- 0.0110 

-0.0035 


0.2771 
0.2297 
0.1356 
0 . 0486 
-0.0046 
-0.0236 
- 0.0211 
-0.0108 
-0.0017 
0.0031 


0.3601 
0 .2604 
0 . 1045 
0 . 0023 
-0.0293 
-0 . 0213 
-0 . 0056 
0 . 0032 
0 . 0044 
0.0023 
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10 

- 0.0138 

- 0.0110 

- 0.0038 

0.0014 

0.0039 

11 

- 0.0108 

- 0.0062 

0.0005 

0.0034 

0.0027 

12 

- 0.0070 

- 0.0021 

0.0027 

0.0033 

0.0010 

13 

- 0.0034 

0.0008 

0.0032 

0.0022 

- 0.0001 

14 

- 0.0005 

0.0024 

0.0026 

0.0009 

- 0.0005 

15 

0.0014 

0.0028 

0.0016 

0.0000 

- 0.0005 

16 

0.0023 

0.0025 

0.0006 

- 0.0004 

- 0.0003 

17 

0.0025 

0.0018 

- 0.0001 

- 0.0005 

- 0.0001 

18 

0.0022 

0.0010 

- 0.0004 

- 0.0004 

0.0001 

19 

0.0017 

0.0003 

- 0.0005 

- 0.0002 

0.0001 

20 

0.0010 

- 0.0001 

- 0.0004 

0.0000 

0.0000 

21 

0.0005 

- 0.0004 

- 0.0002 

0.0000 


22 

0.0000 

- 0.0004 

- 0.0001 

0.0001 


23 

- 0.0002 

- 0.0004 

0.0000 

0.0001 


24 

- 0.0004 

- 0.0003 

0.0001 



25 

- 0.0004 

- 0.0002 

0.0001 



26 

- 0.0003 

0.0000 

0.0001 



27 

- 0.0002 

0.0000 

... 



28 

- 0.0001 

0.0001 




29 

- 0.0001 

0.0001 




30 

0.0000 

0.0001 




31 

0.0000 





32 

0.0001 





33 

0.0001 





34 

0.0001 





35 

• • • 


... 

... 



q n Coefficients 

Auxiliary Data Formula : 

u n+ j = < 5 l u n+ j -1 + ( ? 2 u n + j -2 + ^“ n + j-S + •**' 


n 

6 = 0.01 

6 = 0.02 

6 = 0.05 

6 = 0.10 

6 = 0 . 25 


1 

0 • 9155 

1.0237 

1.1847 

1.3209 

1.5204 


2 

0.4491 

0.4337 

0.3812 

0.3082 

0.1523 

- 

3 

0. 1321 

0.0626 

- 0.0600 

- 0.1750 

- 0.3468 

- 

4 

- 0.0563 

- 0.1293 

- 0.2297 

- 0.2965 

- 0.3471 

- 

5 

- 0. 1442 

- 0.1917 

- 0.2309 

- 0.2294 

- 0. 1700 


6 

- 0. 1619 

- 0.1745 

- 0.1542 

- 0.1060 

- 0.0075 


7 

- 0.1374 

- 0.1199 

- 0.0633 

- 0.0027 

0.0731 


8 

- 0.0937 

- 0.0579 

0.0069 

0.0531 

0.0810 


9 

- 0.0474 

- 0.0065 

0.0452 

0.0654 

0.0523 

- 

10 

- 0.0084 

0 . 0268 

0.0548 

0.0512 

0.0189 

- 

11 

0.0185 

0.0418 

0.0460 

0.0281 

- 0.0031 

- 

12 

0.0330 

0.0425 

0.0293 

0.0078 

- 0.0112 

- 

13 

0.0369 

0 . 0345 

0.0127 

- 0.0045 

- 0.0100 


14 

0.0334 

0 . 0229 

0.0005 

- 0.0091 

- 0.0053 


15 

0.0257 

0.0114 

- 0.0060 

- 0.0083 

- 0.0011 


16 

0.0167 

0 . 0024 

- 0.0079 

- 0.0052 

0.0012 


17 

0.0083 

- 0 . 0034 

- 0.0067 

- 0.0020 

0.0017 

- 

18 

0.0017 

- 0 . 0060 

- 0.0043 

0.0003 

0.0012 

- 

19 

- 0.0027 

- 0.0063 

- 0.0018 

0.0013 

0.0005 

- 


0.0003 

0.0006 

0.0005 

0.0002 

0.0000 

0.0001 

0.0000 


( j > 0 ) 


6 = 1.00 


1.8615 
0.2545 
0.5842 
0.2660 
0.0302 
0.1257 
0.0899 
0.0276 
0.0087 
0.0156 
0.0088 
0.0013 
0.0020 
0.0019 
0.0008 
0.0000 
0.0003 
0.0002 
0 . 0001 
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20 

-0.0049 

-0.0051 

21 

-0.0055 

-0.0034 

22 

-0.0049 

-0.0017 

23 

-0.0037 

-0.0003 

24 

-0.0023 

0.0006 

25 

-0.0010 

0.0010 

26 

-0.0001 

0.0010 

27 

0.0005 

0.0008 

28 

0.0008 

0.0005 

29 

0.0009 

0.0003 

30 

0.0007 

0.0000 

31 

0.0005 

-0.0001 

32 

0.0003 

-0.0002 

33 

0.0001 

-0.0002 

34 

0.0000 

-0.0001 

35 

-0.0001 

-0.0001 

36 

-0.0001 

... 

37 

-0.0001 

... 

38 

-0.0001 

, , , 

39 

-0.0001 

• • • 

40 

• • • 

... 


0.0000 

0.0014 

0.0000 

0.0010 

0.0010 

-0.0002 

0.0012 

0.0005 

-0.0002 

0.0010 

0.0001 

-0.0001 

0.0007 

-0.0012 

• • • 

0.0003 

-0.0002 

• • • 

0.0000 

-0.0002 

• • ♦ 

-0.0001 

-0.0002 

-0.0002 

-0.0001 

-0.0001 

• • • 

• • • 

• m • 

• • • 


6.0 PROTOTYPE ALGORITHM RESULTS 

This section provides an overview of the results obtained from a 
prototype UNIX FORTRAN F77 (i.e., MICROWAY, Inc., NDP FORTRAN-386 
running under MICROPORT'S - SYSTEM V/386 Unix Operating System) . 
This program implements the probabilistic signal reconstruction 
formulations outlined in Sections 3.0 and 5.0. Appendix A, con- 
tains the main statistical nromal probability filtering (SNPF) 
program and supporting normal error probability statistical 
(NEPS.F) subroutine implementations used to produce these re- 
sults. This prototype software is currently configured to con- 
sider only the one minute (coherently integrated) outputs from 
the (software) TYCHO Technology Inc. System. Hence no concensus 
averaging over time is possible under this software implementa- 
tion . 


6.1 GROUND CLUTTER ASSESSMENT 

The first increment of data results consider the first twenty- 
five range gates for beam 1, extracted from the archived data 
file 14130612 . 90P. These lower gate echelons were selected to 
analyze the data filtering effects resulting from local ground 
clutter and non-uniform ground heating mechanisms. Since each of 
these lower range gates have both their raw and filtered dominant 
power spectral peaks near the central frequency, the emphasis for 
the data comparison was shifted to a collation of second largest 
peaks. To aid in the location of these secondary peaks, a trian- 
gle is placed at the filtered peak location in Figures 6.1-1 
through 6.1-24. A summary of the results contained in these Fig- 
ures is provided in Table 6-1. 
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FIGURE 6.1-1 RANGE GATE 2 of BEAM i DEVELOPED FROM 


RANGE GATES 1-3 


Raw Vs Filtered Range Gate Power 

For Range Gate 2 (with 1-3) 



Filter Raw 
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FIGURE 6.1-2 RANGE GATE 3 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 2-4 


Raw Vs Filtered Range Gate Power 

For Range Gate 3 (with 2- 4 ) 



Filter Raw 
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Relative Power 


FIGURE 6.1-3 RANGE GATE 4 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 3-5 


Raw Vs Filtered Range Gate Power 

For Range Gate 4 (with 3- 5 ) 



Power Spectral Frequency 


Filter Raw 
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FIGURE 6.1-4 RANGE GATE 5 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 4-6 


Raw Vs Filtered Range Gate Power 

For Range Gate 5 (with 4- 6 ) 



Filter Raw 
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FIGURE 6.1-5 RANGE GATE 6 OF BEAM i DEVELOPED FROM 
RANGE GATES 5-7 


Raw Vs Filtered Range Gate Power 

For Range Gate 6 (with 5- 7 ) 
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Relative Power 


FIGURE 6.1-6 RANGE GATE 7 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 6-8 


Raw Vs Filtered Range Gate Power 

For Range Gate 7 (with 6- 8 ) 



Filter Raw 
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FIGURE 6.1-7 


RANGE GATE 8 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 7-9 


Raw Vs Filtered Range Gate Power 

For Range Gate 8 (with 7- 9 ) 
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Relative Power 


FIGURE 6.1-8 


RANGE GATE 9 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 8-10 


Raw Vs Filtered Range Gate Power 

For Range Gate 9 (with 8-10) 



Filter Raw 
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FIGURE 6.1-9 


RANGE GATE 10 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 9-11 



Power Spectral Frequency 


Filter Raw 
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Relative Power 


FIGURE 6.1-10 RANGE GATE li OF BEAM 1 DEVELOPED FROM 
RANGE GATES 10-12 



Power Spectral Frequency 


Filter Raw 
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FIGURE 6.1-11 RANGE GATE 12 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 11-13 


Raw Vs Filtered Range Gate Power 

For Range Gate 1 2 (with 1 1 - 13) 



Filter Raw 
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Relative Power 


FIGURE 6.1-12 RANGE GATE 13 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 12-14 
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Relative Power 


FIGURE 6.1-13 RANGE GATE 14 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 13-15 



Power Spectral Frequency 




FIGURE 6.1-14 RANGE GATE 15 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 14-16 




FIGURE 6.1-15 RANGE GATE 16 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 15-17 
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FIGURE 6.1-16 RANGE GATE 17 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 16-18 


Raw Vs Filtered Range Gate Power 

For Range Gate 1 7 (with 16-18) 




FIGURE 6.1-17 RANGE GATE 18 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 17-19 
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FIGURE 6.1-18 


RANGE GATE 19 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 18-20 


Raw Vs Filtered Range Gate Power 

For Range Gate 1 9 (with 18-20) 
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Relative Power 


FIGURE 6.1-19 RANGE GATE 20 OF BEAM 1 DEVELOPED FROM 

RANGE GATES 19-21 
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FIGURE 6.1-20 RANGE GATE 21 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 20-22 


Raw Vs Filtered Range Gate Power 

For Range Gate 21 (with 20- 22 ) 
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FIGURE 6.1-21 RANGE GATE 22 OF BEAM 1 DEVELOPED FROM 

RANGE GATES 21-23 


Raw Vs Filtered Range Gate Power 

For Range Gate 22 (with 21 - 23 ) 



Filter Raw 
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FIGURE 6.1-22 RANGE GATE 23 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 22-24 


Raw Vs Filtered Range Gate Power 

For Range Gate 23 (with 22- 24 ) 
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Relative Power 


FIGURE 6.1-23 RANGE GATE 24 OF BEAM 1 DEVELOPED FROM 
RANGE GATES 23-25 
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FIGURE 6.1-24 RANGE GATE 25 OF BEAM 1 DEVELOPED FROM 


RANGE GATES 24-26 



Power Spectral Frequency 


Filter Raw 
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Given the variability in frequency bin locations and raw data 
variance of Table 6-1, between range gates 9 through 19, a second 
increment of data was computed. This process used two range 
gates on each side of the desired range gate. These results are 
provided in Table 6-2. 


TABLE 6-1 

SUMMARIZED BEAM 1, LOWER RANGE GATE 
RESULTS WITH ONE SURROUNDING RANGE GATE 


RANGE 

GATE 

NUMBER 


FILTERED SECOND LARGEST PEAK 
RELATIVE FREQUENCY 

POWER BIN 

AMPLITUDE 


RAW 

DATA 

VARIANCE 


FILTERED 

DATA 

VARIANCE 


2 

623 . 59 

234 

37.83 

8 .83 

3 

370.86 

144 

30.92 

7 . 13 

4 

424.24 

144 

29.76 

6.69 

5 

255.89 

145 

28.93 

6.42 

6 

150.28 

144 

28.68 

6.29 

7 

67.96 

143 

24.73 

9.79 

8 

47.16 

143 

25.79 

5.67 

9 

24 . 37 

143 

24.39 

5.37 

10 

17.54 

143 

26.28 

5.59 

11 

13 . 38 

117 

26.18 

5.29 

12 

3.48 

116 

24.28 

5.30 

13 

2 . 03 

116 

20.90 

6.96 

14 

1.84 

111 

— 


15 

2 . 07 

161 

23.54 

5.43 

16 

1.93 

89 

22.20 

5.32 

17 

2 .26 

190 

24.50 

5.39 

18 

3 . 02 

89 

21.87 

5.23 

19 

9.26 

89 

18.90 

8 . 29 

20 

110.95 

86 

22.61 

5 .36 

21 

272 . 69 

86 

25.08 

10 .24 

22 

264 . 11 

86 

24.32 

9.68 

23 

141.43 

85 

26.24 

5.56 

24 

56.99 

146 

22.01 

4 . 89 

25 

84 .49 

109 

21.50 

5.28 
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TABLE 6-2 

RESULTS FROM BEAM 1 WITH 
TWO SURROUNDING RANGE GATES 


RANGE 

GATE 

NUMBER 


FILTERED LARGEST PEAK RAW 
RELATIVE FREQUENCY DATA 

POWER BIN VARIANCE 

AMPLITUDE 


FILTERED 

DATA 

VARIANCE 


PROBABLE 

MEAN 

DATA 

ERROR 


9 

2265.29 

129 

10 

1330.42 

128 

11 

894 .85 

128 

12 

695.73 

128 

13 

972.96 

127 

14 

958.98 

126 

15 

1019.04 

126 

16 

1304.51 

126 

17 

3983 . 87 

125 

18 

7537 .51 

125 

19 

- 

125 


26.62 

5.46 

. 0032 

23.05 

5.12 

. 0030 

25.23 

5.31 

. 0033 

22.68 

5.11 

. 0029 

23.70 

4.88 

. 0029 

22.03 

4.705 

. 0026 

22.63 

3 . 14 

. 0030 

22.12 

3.11 

. 0029 

21.15 

3.13 

. 0025 

21.68 

5.05 

. 0029 

20.68 

5.18 

. 0027 


6.2 MAXIMUM PEAK COMPENSATION 


Power spectral estimates were generated for range gates 2 through 
25 by the sequential application of the unbiased weighted mean of 
minimum variance estimator and normal probability filter. The 
question of maximum peak compensation must be considered. This 
requirement is mandated for two distinct reasons. The first 
being that maximum power spectral peaks are not generally symmet- 
ric about their maximum values. Secondly, the normal probability 
filter imposes the constraint that isolated raw data peaks are a 
low probability event. Consequently, these erroneous peaks are 
removed, which tends to aggravate the degree of symmetric peak 
uniformity. 

To compensate for this lack of maximum peak symmetry, Green's 
Theorem can be efficiently used to compute the peak area and 
center of gravity by forming a closed region around the peak from 
the localized piece— wise values. For example, consider the 
illustrative peak configuration of Figure 6.2-1. 


To compute the area of this closed region from Green's Theorem, 

dP(x,y) dQ(x,y) 


[P(x,y)dx + Q ( x , y ) dy ] = 


[ 


dx 


dy 


] dxdy . 


Where, R is the piece-wise closed region formed by traversing the 
indicated points in the counter-clockwise direction and C is the 
simple closed curve of the region boundary formed by the connect- 
ing straight lines, 

L^ = y = m^x + b^ for i = 1 , 2 , ...7. 
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FIGURE 6.2-1 CLOSED REGION AROUND MAXIMUM POWER PEAK 


Maximum Power Peak 




Here, P(x,y) = 0 and Q(x,y) = x so that 


1 

* r 

m^xdx = 


* 

* 


dxdy = Area. 


For the i th line segment, 


A i " 


m^x 


c il 


Yi+l-Vi . , 2 

) ( x i+l 


Xi 2 ) 


x i+l" x i 


Upon summing these segments over i, the total area, A, of the 
closed region becomes: 

1 7 

A = £ (Yi+l-^i) ( x i+l +x i> • 

2 i=l 


It is interesting to note that the center of gravity, E , of 
this area can be determined from Green's Theorem by setting 
P(x,y) = 0 and Q(x,y) = x 2 . 


Therefore , 


m^x 2 dx = 


2xdxdy = 


so that 3!^ has the form 


*i = 


IUj_X 


3 i x i+l 


1 Yi+l’Yi w 3 3 . 

( ) ( x i+l " x i ) 

6 x i+l~ x i 


or 


1 7 


E = 


6A i=l 


2 (Yi+l-Yi) ( x i+l 2 +x i x i+l +x i 2 ) * 


Here, E , is the frequency location of the Doppler Shift asso- 
ciated with the maximum power spectral peak of the selected range 
gate. Figure 6-2-2 contains a graph of the center of gravity 
compensation relative to the maximum range gate power peaks for 
beam 1 associated with range gates 2 through 25. 
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Range Gate 


Frequency Location 
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FIGURE 6.2-2 BEAM Is MAXIMUM PEAK CENTER OF GRAVITY COMPENSATION 





To illustrate the comparative relationship between these computed 
values and measured balloon velocities, the computed and measured 
values are expressed in terms of their distance from the central 
frequency. Figure 6 — 2 — 3 provides a graphic of their central 
frequency comparison. When viewing this comparison, it must be 
remembered that the measured balloon velocities have both a time 
interval and an inertia bias. 


6.3 LINEAR CORRELATION ASSUMPTION OF ADJACENT RANGE GATE 

To assess the statistical assumption of linearly correlated wind 
velocity with respect to adjacent range gates, a series of beam 3 
range gate results were analyzed. Here, beam 3 was selected 
because these power spectral estimates have a much larger veloci- 
ty resolution capability, i.e., 6.7m/s. 

The summarized results of this analysis can be illustrated by 
considering range gate 70 with three surrounding range gates on 
either side of range gate 70. Figure 6.3-1 denotes the raw vs. 
filtered range gate power for this example. Here, the isolated 
raw central data peak is located at frequency bin 130, with a 
power amplitude of 5.09. However, this peak is washed out (can- 
celled) by the lack of non-linear correlations in the spectral 
power amplitudes of the adjacent range gates shown in Figure 6.3- 
2. These non-linear spectral power amplitude correlations are 
presented in Figures 6.3-3, 6.3-4, and 6.3-5 for three, two and 
one surrounding range gates; respectively. Further amplification 
is provided by the computation of the maximum power peak centroid 
values listed in Table 6.3-1. 


TABLE 6.3-1 

COMPARATIVE PEAK CENTROIDS AS A FUNCTION 
OF SURROUNDING RANGE GATES 

RANGE GATE 70 
BEAM 3 


SURROUNDING 

RANGE 

GATES 


COMPUTED 
MAXIMUM PEAK 
CENTROID 


1 

2 

3 


123.91858 

124.50734 

124.72048 


7.0 CONCLUSIONS AND RECOMMENDATIONS 

The specific objectives of this study were: 1) to recommend a 
preferred algorithm for improving atmospheric signals that are 
occasionally obscured by natural and man-made interference, and 
2) to demonstrate that the preferred statistical Normal Probabil- 
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Range Gate 



Measured From Central Frequency 





FIGURE 6.3-1 RAW vs FILTERED RANGE GATE POWER FOR BEAM 1: 
WITH 3 SURROUNDING RANGE GATES 


Raw Vs Filtered Range Gate Power 

For Range Gate 70 (with 67- 73 ) 
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Relative Power 


FIGURE 6 . 3-2 ENLARGEMENT OF FIGURE 6 . 3-1 DATA 


Raw Vs Filtered Range Gate Power 

For Range Gate 70 (with 67- 73 ) 
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Raw Range Gate Data For: 

Range Gates 67, 70, And 73 




Gate 68 — Gate 70 Gate 72 
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< Raw Range Gate Data For: 

Range Gates 68, 70, And 72 




Gate 69 — Gate 70 Gate 71 


Relative Power 



• Raw Range Gate Data For 

Range Gates 69, 70, And 71 




ity Filtering (SNPF) algorithm used in conjunction with the 
weighted mean of minimum variance estimator provides improved 
reliability of wind velocities. 

Based on the analytical data results obtained from a prototype 
computer program for beam 1 and range gates 2 through 25 of the 
archived data file 14130612 . 90P, the computed radical wind veloc- 
ities appear to be in correspondence with similar velocities 
obtained from balloon measurements. See Figure 6-2-3 for a com- 
parison of these results. When viewing this comparison, one 
should realize that the measured balloon velocities suffer from 
inertia bias and wind dynamic effects. 

The recommendations for further study are: 

(1) The present SNPF algorithm has been tested against one 
limited set of wind predictions. Clearly, an extensive 
comparison of the SNPF predictions with all range gates 
and for numerous noise and interference conditions is 
required. 

(2) In addition, a study is required that includes the 
identification of capabilities and limitations of the 
basic wind profiler and balloon measurements, such as 
balloon dynamics in wind. This information will allow 
an understanding of the intricacies of both wind 
measurement methods and potential explanations of the 
differences in the predictions. 
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APPENDIX A 


PROTOTYPE FORTRAN SOURCE 
CODE 


(Statistical Normal Probability Filtering (SNPF) Programs) 


A-0 



10 


15 


20 


25 


program neps_filter 
implicit real*8 (a-h,o-z) 

dimension jmode(3) , araw(258 , 11) ,var(ll) ,a(ll) ,xbar(256) 
character*28 str,st 
character*13 f ilel , f ile2 , f ile3 
character*2 cr 
character* 1 yn,ny 
data jmode/0, 4816 , 9632/ 
cr = i mi 
vtol = 1 . 5d-10 
iup = 112 
ny = * y 1 

I***********************************' 

** Enter Name Of Input Data File **' 




’)-! 


Name Is: ',filel 

status= ' old ' , f orm= ' formatted ' , 


print* , 
print* , 
print* , 
print* , 

read(*,900) file2 
1 = index ( file2 , 1 
filel = file2(:l) 
if(i.gt.l2) then 
print* , ' 1 

print* ,' Named Input Data File Must Have 12 Or Less' 
print*, 'Alphanumeric Characters Including The Extension.' 
print*, 'Do You Wish To Enter Another File Name (y,n)?' 
read(*,901) yn 
if(yn.eq.ny) goto 5 
goto 999 
endif 

print* ,' Requested Input Data File 
open (8 , f ile=f ilel , access= ' direct ' 
lerr = 35, iostat=if lag, recl=79) 
print* , ' ' 

print* ,' Select Beam Mode (1,2,3)?' 
read* , mode 

print* ,' Selected Beam Mode Equals: ',mode 

if ( (mode . ge. 1) . and. (mode. le . 3 ) ) then 
print* , ' ' 

print* ,' Enter Desired Range Gate (1-112).' 
read*,j 

irange = 2000+ ( j — 1 ) *150 

print*, 'Requested Range Gate Equals: ',j 

print*, 'At Height: ', irange,' (Meters)' 

if ( ( j . It. 1) .or. ( j .gt. 112) ) then 

print* ,' Entered Range Gate Out of Bound.' 
print*, 'Must Be In The Closed Interval [1,112]' 
print*, 'Do You Wish To Renter This Number (y,n)?' 
read*,yn 

if(yn.eq.ny) goto 15 
goto 999 
endif 
print* , ' ' 

print* ,' Enter Number Of Surrounding Range Gates (1-5)' 
print*, 'To Be Considered In The Statistical Data Filter' 
print* ,' Smoothing Process.' 
read*,i 

print*, 'Selected Number Of Surrounding Range Gates Equals: 
if ( (i . It . 1) . or . (i.gt . 5) ) then 

print*, 'Selected Number Of Range Gates Out Of Bound.' 
print*, 'Must Be In The Closed Interval [1,5]' 
print*, 'Do You Wish To Reenter This Number (y,n)?' 
read* , yn 

if(yn.eq.ny) goto 20 
goto 999 
endif 

else „ , , 

print* ,' Selected Beam Mode Is Out Of Bound. 

A- 1 


print*,' Must Be In The Closed Interval [1,3]' 
print*, 'Do You Wish TO Reenter This Number (y,n)?' 
read* , yn 

if(yn.eq.ny) goto 10 
goto 999 
endif 
il = j-i 
ih = j + i _ 
iouth = i 

if (ih.gt. iup) iouth = iup-j 
ioutl = i 

if(il.lt.l) ioutl = j-1 
k = minO (ioutl, iouth) 
if(k.lt.i) then 

print*, 'Requested Number Of Surrounding Range' 
print* , 'Gates Cannot Be Supported Using The' 
print*, 'Currently Specified Range Gate Value. ' 
print*, 'This Number Must Be Less Than Or Equal' 
print* , ' To ' , k 
goto 25 
endif 
pir i nt * / * * 

print * ' Screen Print Of Incrementally Processed Data (y,n)?' 
read* , yn 
iwrt = 0 

if(yn.eq.ny) iwrt = 1 
igate = 1 
kk = 256 

do istart = il,ih 
print* , ' ' 

if ( iwrt . eq. 1 ) print*, 'Raw Data For Range Gate ', istart 
icount = 1 

iend = 43*istart+jmode (mode) 
ibegin = iend-42 
do 1 = ibegin, iend 

read (8 , 902 , rec=l , iostat=if lag, err=30) (a(i) , i=l, 6) 

do i = 1 , 6 

araw( icount, igate) = a(i) 
if ( ( iwrt . eq. 1) .and. ( icount . le . kk) ) then 
print*,' ', icount,' ',a(i) 

endif 

icount = icount+1 
end do 
end do 

igate = igate+1 
end do 

igate = igate/2 

c utilize NEPS to obtain raw power spectral range gate 

c data error variances 

k = 1 
kk = 256 
do i = il,ih 
do j = 1 , kk 

c temporarily load the raw power spectral range gate 

c data into the xbar array, 

xbar ( j ) = araw ( j ,k) 
end do 

call neps (xbar, kk, segma , pme , ierr ) 
var(k) = segma*segma 

print*, 'Data Variance/Probable Error For Range Gate ',i 
print*,' Variance = ' ,var(k) 

print*,' Probable Error = ',pme 

k = k+1 

end do . . 

compute unbiased weighted mean of minimum variance weighting 
coefficients and variance estimates. 
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ngates = ih-il+1 
sum = 0 . OdO 
do j = 1, ngates 

if (var(j) .le.vtol) var(j) = vtol 
sum = sum+ ( 1 . OdO/var ( j ) ) 
end do 

if (iwrt . eq. 1) then 
print* , ' ' 

print*, 'Unbiased Range Gate Coefficient Weights:' 
endif 
ip = il 
vxbar =0.0d0 
do j = 1, ngates 

a(j) = 1. OdO/ (var ( j ) *sum) 
vxbar = vxbar+a ( j ) *a ( j ) *var ( j ) 
if ( iwrt . eq. 1) then 

print*, ip, 'th. Range Gate Coefficient = ',a(]) 
ip = ip+1 
endif 

end do t . 

c compute weighted mean of minimum variance estimates 

do j = 1 , kk 

xbar(j) = O.OdO 
do i = 1, ngates 

xbar(j) = xbar(j)+a(i)*araw(j,i) 
end do 
end do 

i = ( il+ih) /2 

print* , ' ' . , 

print*, 'Weighted Mean Of Minimum Variance Equals: ', vxbar 

print*, ' Iterated Data Error Variance: (range gate ',i,| )' 

c statistically filter the weighted mean of minimum variance 

c estimates down to a data error variance of vxbar 

call neps (xbar, kk, segma, pme, ierr) 
vfmmv = segma*segma 
print*,' ', vfmmv 

do while (vfmmv. gt. vxbar) 

call neps ( xbar, kk, segma, pme, ierr) 
vfmmv = segma*segma 
print*,' ', vfmmv 

end do 
print* , ' 1 

c compute frequency index at maximum power 

xmax = xbar(l) 
jmax - 1 
do j = 2 , kk 

if (xbar(j) .gt.xmax) then 
xmax = xbar(j) 
jmax = j 
endif 
end do 

dopsmp = 128-jmax 

print* ,' Doppler Shift At Maximum Power Peak = ', dopsmp 
print* ,' Power At Maximum Peak = ' , xmax 
print* , ' ' 

c compute center of gravity of maximum power peak using 

c green's theorem, 

ip = jmax+2 
a(l) = dble(ip) 
var(l) = O.OdO 
do j = 2,6 

a ( j ) = dble(ip) 
var(j) = xbar (ip) 
ip = ip-1 
end do 
a ( 7 ) = a ( 6) 
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var(7) = O.OdO 
area = 0 . OdO 
vel = area 
do j = 1 , 6 

j 2 = ^ +1 

temp = var ( 32 ) -var ( 3 ) 

area = area+temp* (a ( j 2 ) +a ( j ) ) . . 

vel = vel+temp*(a(j 2 ) *a ( j 2 ) +a ( j ) *a (j 2 ) +a (j ) *a (j ) ) 
end do 

area = 0.5d0*area 
vel = ( 1 . OdO/ ( 6 . OdO* area) ) *vel 
print*, 'Radial Wind Velocity = ',vel 
print* , ' ' 

if (iwrt . eq. 1 ) then , . , 

print* Filtered Mean Of Minimum Variance Estimators 

do j = 1 , kk 

print*, ' ' , j , ' ' , xbar ( j ) 

end do 

c write ASCII comma and "" delimited file for QUATTO PRO plotting. 

1 = index ( filel ) 
f ile 2 = filel ( : 1 ) // ' pit 1 
f ile 3 = filel ( : 1) // 1 tern' 
print*, ' ' 

print* ,' Currently Generating File For' 
print*, 'Later QUATTO PRO Plotting.; 
print*, 'Plot File Name Equals: ' ,file2 

open ( 9 , f ile=f ile3 , access= ' direct ' , status= ' scratch ' , 
lform= ' formatted ', recl=26) 

open ( 10 , f ile=f ile 2 , access= ' sequential ' , status= unknown , 
lform= ' unformatted ' ) 

c write required QUATTO PRO plotting series, 

c axes and text legends 

i = (il+ih )/2 
st ( 1 : 2 ) = cr 
write ( 10 , 906) st ( 1 : 2 ) 
write ( 10 , 907 ) i, il , ih, st ( 1 : 2 ) 
write (10,908) st(l:2) 

write(10, 909) st(l:2) 

write(10, 910) st(l: 2 ) 

write ( 10 ,-911) st(l:2) 

do j = 1 , kk 

write (9 , 903 , rec=j ) j ,xbar(j) ,araw(] ,igate) 
end do 
do j = 1 , kk 

read (9 , 904 , rec=j ) st(:26) 

str ( : 28 ) = st(:26)//cr 
write (10,905) str(:28) 

end do 



close (unit=10) 
close (unit=8 ) 
goto 999 
print* , 1 ' 

if ( if lag. eq. -1) then 
print* , ' ' 

print*, 'End-Of-File Condition 
print*, 'On Input Data File: ' 

goto 999 


Has Ocurred' 
filel 


else . 

print*, 'Read Operation on Input Data File:' 
print*,' (', filel,')' 

print*, 'Has Been Terminated With A Standard 
print* ,' FORTRAN Run-time Error Of: ' , if lag 

goto 999 

print*, 'Unknown Problem With Requested Input Data File. ' 
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print*,' ('jfilel,')' 

print*, 'Must Terminate File Processing.' 

900 format (al3) 

901 format (al) 

902 format (6 (el3 . 6) , lx) 

903 format ( i8 , 2 ( 1 , ' , f 8 . 2 ) ) 

904 format (a26) 

905 format (a28) 

906 format (34h"Raw Vs Filtered Range Gate Power", a2) 

907 format (16h"For Range Gate ,i3,7h (with , i3 , lh-, 13 , 3h )",a2) 

908 format (26h"Power Spectral Frequency ", a2 ) 

909 format (16h"Relative Power",a2) 

910 format (8h"Filter" , a2 ) 

911 format (5h"Raw",a2) 

999 stop 

end 
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SUBROUTINE NEPS (Y,N, SEGMA , PME , IER) 

IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION A (1072) , Y(1000) ,BU(39) ,AJ(173) ,JA(18) ,AK(153) f KA(12) 
PARAMETER (FACT = 0.6744385D0) 

NEPS LAST MODIFIED 12-17-91 

DATA KA/1, 35, 36, 66, 67 , 93 , 94 , 117 , 118 , 137 , 138, 153/ 
############################################################## 12 

DATA JA/100, 1,39, 50, 40, 74, 20, 75, 104, 10, 105, 131, 4, 13 2, 154, 1, 155, 

$ 173/ 

##############################################################18 
DATA AK/0 . 157 0D0, 0 . 14 82 DO, 0. 12 54 DO , 0 . 0948D0 , 0 . 062 8 DO , 0 . 034 IDO , 

$ 0. 0117D0, -0 . 0035D0, -0. 0119D0, -0. 0148D0, -0. 0138D0, -0 . 0108 DO, 

$ -0. 0070DO, -0 . 003 4 DO , -0 . 0005D0, 0. 0014D0, 0 . 0023 DO, 0 . 0025D0 , 

$ 0 . 002 2D0 , 0. 0017D0 , 0 . 0010D0, 0 . 0005D0 , 0. 0000D0 , -0 . 000 2 DO , 

$ -0. 0004 DO, -0. 0004 DO, -0.0003 DO, -0.0002 DO, -0.0001D0, -0.0001 DO, 

$ 0 . OOOODO , 0 . OOOODO, 0 . 000 IDO , 0 . 000 IDO , 0 . 0001D0, 

##########################################35################### 

$ 0 . 17 6 9 DO , 0 . 164 4 DO , 0 . 1329D0 , 0 . 0928D0 , 0 . 053 6D0 ,0.0216D0, 

$ -0 . 0004 DO, -0 . 0125D0 , -0 . 0165D0, -0 . 0151D0 , -0 . 0110D0 , -0 . 006 2 DO , 

$ -0 . 002 IDO , 0. 0008 DO, 0.0024 DO, 0.0028D0 , 0. 0025D0, 0. 00 18 DO , 

$ 0 . 00 10 DO ,0.003D0,-0. 000 IDO , -0 . 0004 DO , — 0 . 00 04 DO , -0 . 00 04 DO , 

$ -0. 0003D0, -0.0002D0,0. OOOODO, 0. OOOODO, 0 . 0001D0 , 0 . 0001D0 , 

fttSaSeeJ?####################################################### 

$ 0.2076D0, 0. 1873D0, 0. 1397D0 , 0 . 084 2D0 , 0 . 0359D0 , 0 . 002 7 DO , 

$ -0. 0145D0 , -0. 019 IDO, -0 . 016 IDO, -0 .0099D0,-0.0038D0,0.0005D0, 

$ 0. 0027 DO , 0 . 00 3 2 DO , 0. 0026D0 , 0 . 0016D0 , 0 . 0006D0, -0 . 000 IDO , 

$ -0. 0004 DO, -0. 0005 DO, -0.0004 DO, -0.0002D0, -0.000 IDO, 0.0000 DO, 

#######################93####################################### 

$ 0 . 2 3 47 DO 0 . 2056D0 , 0 . 14 12 DO ,0.0721D0,0. 019 IDO ,-0.0110D0, 

$ -0. 0211D0, -0 . 018 6D0 , -0 . 0110D0 , -0 . 003 5D0 , 0 . 00 14 DO , 0 . 00 3 4 DO , 

$ 0. 0033D0, 0. 0022D0, 0. 0009D0, 0. OOOODO, -0. 0004D0, -0. 0005D0 , 

$ -0. 0004 DO, -0 . 0002 DO , 0 . OOOODO , 0 . OOOODO , 0 . 0001D0 , 0 . 000 IDO , 

#####################################################117####### 

$ 0.22771D0, 0.2297D0, 0. 1357D0, 0 . 0486D0 , -0 . 0046D0 , -0 . 0236D0, 

$ 0. 02 11D0, -0. 0108D0, -0.0017D0, 0.003 IDO, 0.003 9D0, 0.002 7D0, 

$ 0. 0010D0, -0. 0001D0, -0.0005 DO, -0. 0005D0, -0.0003 DO, -0.0001D0, 

##############137############################################## 

$ 0 . 3 60 IDO , 0 . 2 6 04 DO , 0 . 1045D0, 0 . 002 3 DO , -0 . 029 3D0 , -0 . 02 13 DO , 

$ -0 . 00 5 6 DO , 0 . 00 3 2 DO , 0 . 004 4 DO , 0 . 002 3 DO , 0 . 00 03 DO ,-0. 0006D0, 

S -0. 0005D0 , -0 . 00 02 DO , 0 . OOOODO , 0 . 000 IDO/ 

################################# 153########################173 
DATA AJ/0.9155D0,0.4491D0,0.1321D0, -0.0563 DO, -0.1442 DO, -0.1619 DO, 

$ -0. 1374 DO, -0. 0937 DO, -0.04 74 DO, -0.0084 DO, 0. 018 5D0 , 0 . 0330D0 , 

$ 0. 0369 DO, 0. 0334 DO, 0. 0257 DO, 0.0167 DO, 0.0083 DO, 0.0017D0, -0.0027 DO, 

$ -0. 004 9 DO, -0. 0055D0, -0 . 0049D0, -0 . 003 7 DO, -0 . 002 3 DO , -0 . 0010D0 , 

$ — 0 . 000 IDO ,0.0005D0,0.0008D0,0. 00 09 DO , 0 . 00 07 DO , 0.0005D0, 0. 0003D0, 

$ 0. 000 IDO, 0. OOOODO, -0 . 0001D0 , -0 . 0001D0 , -0 . 0001D0 , -0.000 IDO, 

$ -0. 0001D0, 1. 02 37D0, 0. 43 37D0, 0. 0626D0, -0. 1293D0, -0. 19 17 DO, 

$ -0.1745 DO, -0.1199 DO, -0.0579 DO, -0. 0065D0 ,0.0268 DO, 0.0418 DO, 

$ 0.0425D0,0.0345D0,0.0229D0,0.0114D0,0. 00 2 4 DO ,-0.0034DO,-0.0060DO, 
$ -0.0063D0,— 0.0051D0, -0 . 003 4D0 , -0 . 0017D0 , -0 . 0003D0 , 0. 0006D0, 

$ 0. 00 10D0, 0. 0010D0, 0. 0008 DO, 0. 0005D0, 0. 0003 DO, 0.0000 DO, -0.0001D0, 

$ -0. 0002 DO, -0. 0002 DO, -0.0001 DO, -0.0001D0, 1.1847 DO, 0.3812 DO, 

$ -0. 0600D0, -0 . 2297D0 , -0 . 2 3 09D0 , -0 . 154 2D0 , -0 . 063 3 DO , 0 . 00 69 DO , 

$ 0.04 52 DO, 0.0548 DO, 0.0460 DO, 0.0293 DO, 0.0127 DO ,0.0005 DO, -0.0060 DO, 

$ -0. 0079D0, -0. 0067 DO, -0. 0043 DO, -0. 00 18D0,0. OOOODO, 0.0010 DO, 

$ 0 0012 DO, 0.0010 DO, 0.0007 DO, 0.0003 DO, 0.0000 DO, -0.0001 DO, -0.0002 DO, 

$ — C).0002D0,-0.0001D0, 1 . 3209D0 , 0 . 3082D0 , — 0 . 1750D0 , -0 . 2 965D0 , 

$ -0.2294 DO, -0. 1060D0, -0.0027 DO, 0. 053 IDO ,0.0654 DO, 0.0512 DO, 

$ 0. 02 8 IDO , 0 . 007 8 DO, -0. 004 5D0 , -0 . 009 IDO , -0 . 008 3 DO , -0 . 0 052 DO , 

$ -0.0020 DO, 0.0003 DO, 0.0013 DO, 0.0014 DO, 0.0010DO,0.0005DO,0.0001DO, 
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$ -0 . 00 02 DO , — 0 . 00 02 DO , — 0 . 000 2 DO , -0 . 000 1D0 , 1 . 5 2 04 DO , 0 . 152 3D0 , 

$ - 0 . 34 68D0,-0. 347 1D0, -0.1700D0,-0. 0075D0, 0.073 1D0, 0.08 10D0, 

$ 0. 052 3 DO, 0.0 189 DO, -0.003 1D0, -0 . 0112D0 , -0 . 0100D0 , -0.0053 DO, 

$ -0 . 0011D0 , 0 . 00 12 DO , 0 . 00 17 DO , 0 . 00 12 DO , 0 . 0005D0 ,0.0000D0, 

$ -0.0002D0, -0. 0002D0, -0. 0001D0, 1 . 8615D0 , -0 . 2545D0 , -0 . 5842D0, 
$ -0. 2660D0, 0 . 03 02 DO , 0 . 1257D0 , 0 . 0899D0 , 0 . 0276D0 , -0 . 00 87 DO , 

$ -0 . 015 6D0 , -0.0088 DO, -0.0013 DO, 0. 002 0D0 , 0 . 0019D0 , 0 . 0008D0 , 

$ 0.0000D0, -0.0003 DO, -0. 0002D0 , -0 . 000 IDO/ 

DETERMINE MINIMUM DATA OBSERVATION VALUE SO THAT 

THE DATA CAN BE TRANSFORMED, ( I . E ., Y '( I) = Ln ( Y ( I ) -YMIN) , 

SO THAT ITS MEASURE OF PRECISION HAS NEARLY THE SAME 
VALUE FOR ALL OBSERVATION VALUES. 

YMIN = Y ( 1 ) 

DO I = 2 , N 

YMIN = 0 . 5D0* (YMIN+Y (I) -DABS (YMIN-Y (I) ) ) 

END DO 

YMIN = YMIN-2 . 0D-08 

CHECK ADEQUACY OF OBSERVATION DATA SAMPLE. 

IF(N.GT. 4) THEN 
SUM = 0.0D0 
DO I = 2 , N-2 

R = (Y(I) -YMIN)/ (Y(I+1) -YMIN) 

S = (Y(I+2) -YMIN)/ (Y(I-l) -YMIN) 

A ( 1+2 ) = DLOG (S) +3 . ODO*DLOG (R) 

SUM = SUM+A ( 1+2 ) 

END DO 

COMPUTE THIRD DIFFERENCE MEAN. 

SUM = SUM/DBLE (N-3 ) 

COMPUTE TRANSFORMED PROBABLE THIRD DIFFERENCE ERROR 
FROM THE NORMAL PROBABILITY INDEX OF PRECISION. 

PE = 0.0D0 
DO I = 4 , N 

R = A (I) -SUM 
PE = PE+R*R 
END DO 

Cl = DBLE (N-3 ) 

C2 = DSQRT(PE) *FACT 
PE = C2/DSQRT (Cl-1 . 0D0) 

COMPUTE MEAN PROBABLE SMOOTHING ERROR REDUCTION. 

PME = C2*PE/ (DSQRT (C1*C1-C1) *C1) 

COMPUTE PROBABLE ERROR STANDARD DEVIATION. 

SEGMA = PE/ FACT 
ELSE 

INSUFFICIENT DATA SAMPLE FOR STATISTICAL DATA SMOOTHING. 
MUST HAVE GREATER THAN FOUR, (4), DATA OBSERVATIONS. 

IER = 1 
RETURN 
ENDIF 
IER = 0 

USING THE FUNDAMENTAL THEOREM OF INDUCTIVE PROBABILITY, 

MAP THE OBSERVED THIRD DIFFERENCE PROBABLE DATA ERRORS 


A-7 



nnn ooo oonno ono oooooo non 


INTO THE SET OF SMOOTHING COEFFICIENTS. 


DO 65 I = 1,16,3 

EPS = 1 . ODO/DBLE ( JA ( I ) ) 

IF ( (PE.LE.EPS) .OR. (I.EQ. 16) ) THEN 

EXTRAPOLATE UPPER AUXILIARY POINT DATA. 

L = 1+1 

M = JA(L+1) -JA(L) +1 
IF(M.GT.N) THEN 

STATISTICAL DATA SAMPLE NOT LARGE ENOUGH TO PROVIDE 
HIGH CONFIDENCE DATA SMOOTHING OR AUXILIARY DATA POINT 
GENERATION. NOTICE, THE MINIMUM SAMPLE SIZE IS RETURNED 
IN IER. 

IER = M 
M = N 
ENDIF 

DETERMINE NUMBER OF REQUIRED AUXILIARY DATA POINTS. 

IP = (2*1+1) /3 

NK = KA(IP+1) -KA(IP) 

DO 10 K = 1 , M 

BU (K) = DLOG ( Y (K) -YMIN) 

LO CONTINUE 

GENERATE AUXILIARY DATA POINTS WHICH SATISFY THE A PRIORI 
PROBABILITY HYPOTHESIS THAT THE SUPPLIED DATA OBSERVATIONS 
HAVE A NORMAL ERROR LAW DISTRIBUTION. 

Ml = M-l 

DO 25 K = NK, 1,-1 
LI = JA ( L) 

SUM = AJ ( LI ) * BU ( 1 ) 

LI = LI— 1 
DO 15 J = 2,M 

SUM = SUM+AJ (Ll+J) *BU (J) 

15 CONTINUE 

A (K) = SUM 

REASSIGN UPPER AUXILIARY POINTS. 

DO 20 LI = Ml, 1, -1 
BU ( Ll+1 ) = BU (LI) 

20 CONTINUE 

BU ( 1 ) = SUM 
25 CONTINUE 

EXTRAPOLATE LOWER AUXILIARY POINT DATA. 

LI = 1 

DO 30 K = N, N-M+l , -1 

BU (LI) = DLOG (Y(K) -YMIN) 

LI = Ll+1 
30 CONTINUE 

NK1 = NK+1 

DO 4 5 K = N+NK1 , 2 *NK+N 
LI = JA ( L) 

SUM = AJ (LI) *BU ( 1) 

LI = Ll-1 
DO 35 J = 2 , M 

SUM = SUM+AJ (Ll+J) *BU(J) 

35 CONTINUE 
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A (K) = SUM 

REASSIGN LOWER AUXILIARY POINTS. 

DO 40 LI = Ml, 1,-1 
BU (Ll+1) = BU(L1) 

0 CONTINUE 

BU ( 1) = SUM 
5 CONTINUE 

LOAD OBSERVATION DATA INTO THE A(J) ARRAY FOR 
J = NK+1 , NK+2 , . . . , NK+N . 

DO 50 K = 1,N 

A (NK+K) = DLOG ( Y (K) -YMIN) 

0 CONTINUE 

IMPOSE CONSERVATION THEOREM SMOOTHING CONDITIONS: 

1 . ) THE SUM OF RAW DATA OBSERVATIONS EQUALS THE 
SUM OF SMOOTHED DATA OBSERVATIONS. 

2. ) THE RAW DATA MOMENTS OF ORDERS ONE, TWO, AND 
THREE EQUAL THE MOMENTS OF THE SMOOTHED DATA. 

3 . ) THE SMOOTHED AND UN-SMOOTHED DATA OBSERVATION 
GRAPHS HAVE THE FOLLOWING CHARACTERISTICS: 

a) Equal areas, 

b) Same center of gravity, and 

c) Identical moments of inertia 
about any line parallel to the 
observation ordinates, e. g. , the 
Y-axis . 

LI = KA(IP) 

DO 60 J = NK1 , NK+N 
SUM = AK (LI) *A ( J) 

DO 55 K = 1 , NK 

SUM = SUM+AK ( Ll+K) * ( A ( J+K) +A ( J— K) ) 

55 CONTINUE 

Y(J-NK) = DEXP ( SUM) +YMIN 
- c Y(J-NK) = SUM 

60 CONTINUE 

RETURN 
ENDIF 

65 CONTINUE 
END 
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